% relation between B_mk_in_analytical_formulas = B_mk_in_this_function./(Qk.*besselj(m, Qk))

function [B_m1 B_m2 B_m1_1 B_m2_1] = func_B_mk_l(q, m, CONSTS)

    k0    = CONSTS.k0;
    a     = CONSTS.a;
    eta   = CONSTS.eta;
    eta_a = CONSTS.eta_a;
    c     = CONSTS.c;
    
    p = func_p(q, CONSTS); 
    [q1 q2] = func_qk_from_p_inv(p, CONSTS);
    [H_m1 H_m2] = func_H_m_l(q, m, CONSTS);
    [J_m1 J_m2 J_m1_1 J_m2_1] = func_J_mk(q, m, CONSTS);
    [n1 n2] = func_nk_from_p(p, CONSTS);
    Q = k0*a*q;
    Q1 = k0*a*q1;
    Q2 = k0*a*q2;
    
    bessel_m_1 = besselj(m, Q1);
    bessel_m_2 = besselj(m, Q2);
    
%     B_m1 = (4*pi/c)*(1./(Q1.*bessel_m_1)) .* (k0*a) .* ...
%            (n2.*(eta/eta_a).*J_m2.*H_m2 + (eta/eta_a)*p.*(m./((Q).^2)).*J_m2_1 - ...
%             n2.*H_m2.*H_m2 + ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n2);
%             
%     B_m2 = (4*pi/c)*(1./(Q2.*bessel_m_2)) .* (k0*a) .* ...
%            (-n1.*(eta/eta_a).*J_m1_1.*H_m2 - (eta/eta_a)*p.*(m./((Q).^2)).*J_m1 + ...
%             n1.*H_m2.*H_m2 - ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n1);
%         
%     B_m1_1 = (4*pi/c)*(1./(Q1.*bessel_m_1)) .* (k0*a) .* ...
%              (n2.*(eta/eta_a).*J_m2.*H_m1 + (eta/eta_a)*p.*(m./((Q).^2)).*J_m2_1 - ...
%              n2.*H_m1.*H_m1 + ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n2);
%         
%     B_m2_1 = (4*pi/c)*(1./(Q2.*bessel_m_2)) .* (k0*a) .* ...
%              (-n1.*(eta/eta_a).*J_m1_1.*H_m1 - (eta/eta_a)*p.*(m./((Q).^2)).*J_m1 + ...
%              n1.*H_m1.*H_m1 - ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n1);

    B_m1 = (4*pi/c)*(k0*a) .* ...
           (n2.*(eta/eta_a).*J_m2.*H_m2 + (eta/eta_a)*p.*(m./((Q).^2)).*J_m2_1 - ...
            n2.*H_m2.*H_m2 + ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n2);
            
    B_m2 = (4*pi/c)*(k0*a) .* ...
           (-n1.*(eta/eta_a).*J_m1_1.*H_m2 - (eta/eta_a)*p.*(m./((Q).^2)).*J_m1 + ...
            n1.*H_m2.*H_m2 - ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n1);
        
    B_m1_1 = (4*pi/c)*(k0*a) .* ...
             (n2.*(eta/eta_a).*J_m2.*H_m1 + (eta/eta_a)*p.*(m./((Q).^2)).*J_m2_1 - ...
             n2.*H_m1.*H_m1 + ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n2);
        
    B_m2_1 = (4*pi/c)*(k0*a) .* ...
             (-n1.*(eta/eta_a).*J_m1_1.*H_m1 - (eta/eta_a)*p.*(m./((Q).^2)).*J_m1 + ...
             n1.*H_m1.*H_m1 - ((p.^2)./eta_a).*((m^2)./((Q).^4)).*n1);
         
 
end 
